function [NewRegionNum, NewRegions, NewRegionSampleId] = Partition_NP_BoxConstraint(Region, NewRegionNum, MaxPartitionDepth, SampleSet)
% the partition action for box constrainted problems (Nested Partition)
% INPUT
%   d: problem dimension
%   Region: 1 -by- (2+2d) vector, [partition depth, partitionable, d1-start, d1_end, d2-start, d2_end, ...]
%       partition depth = 0 for the whole domain
%       partitionable: scalar, 1: partitionable, 0: non-partitionable
%   NewRegionNum: scalar, number of new partitioned regions
%   MaxPartitionDepth: scalar, the maximal partition depth
%   SampleSet: ~ -by- (2+d) matrix, [SampleID, objective value, x1, x2, ...]
% OUTPUT
%   NewRegionNum: scalar, number of new partitioned regions
%   NewRegions: RegionNum -by- (2+2d) matrix, the new partitioned regions
%   NewRegionSampleId: ~ -by- 1 cell, the sample id that belongs to each region

%% PARTITION SUBREGIONS
d = (length(Region)-2)/2; % the problem dimension
x_id = (1:d).*2+1; % the start position of each dimension
NewRegions = repmat(Region,[NewRegionNum,1]);
NewRegions(:,1) = Region(1,1)+1; % the partition depth
if NewRegions(1,1)>=MaxPartitionDepth
    NewRegions(:,2) = 0; % non-partitionable
end
[~, longest_d] = max(Region(1,x_id+1)-Region(1,x_id)); % the longest edge
longest_d_id = x_id(longest_d); % the position of this dimension
longest_d_length = (Region(1,longest_d_id+1)-Region(1,longest_d_id))/NewRegionNum; % the new length of this dimension
NewRegions(:,longest_d_id+1) = Region(1,longest_d_id)+(1:NewRegionNum)'.*longest_d_length;
NewRegions(2:end,longest_d_id) = NewRegions(1:(end-1),longest_d_id+1);

%% MODIFY SAMPLE SET ID
NewRegionSampleId = cell(NewRegionNum,1);
if ~isempty(SampleSet)
    for i = 1:NewRegionNum
        belongs = (SampleSet(:,2+longest_d)>=NewRegions(i,longest_d_id))...
            .*(SampleSet(:,2+longest_d)<NewRegions(i,longest_d_id+1));
        NewRegionSampleId{i} = SampleSet(belongs==1,1);
    end
end

end